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We study the impurity effects on the Caroli-de Gennes-Matricon (CdGM) states, particularly on the level spacings in 
a vortex core in a topological s-wave superconductor (SC) by two methods, numerical and analytical. The topological 
s-wave SC belongs to the same class as a chiral p-wave SC, and thus there are two inequivalent vortices in terms of any 
symmetry operation. We take into account this inequivalence and numerically calculate the scattering rates based on an 
improved version of the Kopnin-Kravtsov (iKK) scheme, which enables us to treat the discrete levels in the presence 
of white-noise disorder. We also construct the Andreev equation for the topological s-wave SC and obtain the Andreev 
bound states analytically. We use a correspondence between the wave functions for the Bogoliubov-de Gennes equation 
and the Andreev equation in the iKK scheme and deduce the formula of scattering rates described by the wave function 
for the Andreev equation. With this formula, we discuss the origin of impurity scattering rates for CdGM states of the 
topological s-wave SC and the dependence on the types of vortex related to the inequivalence. 


1. Introduction 

An implementation of topological quantum computation 
(TQC) with low decoherence is a promising application of 
topological superconductors (TSCs ). 1 ’ 21 In addition to this 
fascinating application, recent intensive studies of TSCs have 
been motivated by fundamental interest in the appearance of 
the Majorana fermion in condensed matter physics . 3 - 41 The 
Majorana fermion exists as a topologically protected zero- 
energy bound state near the edge or around the topologi¬ 
cal defects in TSCs and this zero-energy state obeys a non- 
Abelian statistics. The implementation of TQC is realized via 
the braiding operation among the degenerated Majorana zero- 
energy states. 1,5 ' It is necessary to perform the braiding oper¬ 
ation adiabatically in order to avoid nonadiabatic transitions 
of zero-energy states. The typical operation time should be 
longer than the time given by Ti divided by the level spacing 
between the Majorana state and the first excited state (we call 
this condition the “adiabatic condition”). It is thus important 
that the level spacing, i.e., the energy of the first excited state, 
is stable against disorder as well as the zero-energy state also 
in terms of the signal strength of measurements. 6 - 7 ' In this 
work, we focus on the Majorana state in a vortex core and in¬ 
vestigate the impurity effects on the bound states in a vortex 
core called Caroli-de Gennes-Matricon (CdGM) states. 8 ' 

A chiral p-wave SC 1,9 ~ U ' is a typical two-dimensional TSC 
that belongs to the Bogoliubov-de Gennes (BdG) class with 
broken time-reversal symmetry, called “class D”. This class 
is one of ten symmetry classes obtained by Altland and Zirn- 
bauer on the basis of random matrix theory, 12 ' and all possible 
TSCs have been classified into these classes by Schnyder et 
al . 13-141 if we do not take into account the crystal symmetry. 
A TSC in class D has two inequivalent single vortices in terms 


of symmetry operations: the inequivalence is characterized by 
the relative sign of chirality and vorticity. There is a remark¬ 
able difference between the two inequivalent vortices, accord¬ 
ing to an earlier work concerning impurity effects on low- 
energy states localized in a vortex core in a chiral p-wave SC 
by means of quasiclassical theory. 15 - 16 ' The scattering rates 
of the localized states are almost the same as those in the 
normal states when the relative sign is positive. On the other 
hand, the scattering rates are vanishingly small when the sign 
is negative. This is a consequence of coherence effects and ro¬ 
tational symmetry. However, Sr 2 Ru 04 , a candidate for a chi¬ 
ral p-wave SC, has a point group symmetry C 4 in the crys¬ 
tal structure 171 and we could not expect to observe the effect 
of rotational symmetry. Other possible candidates for class D 
TSCs are, for example, the v = 5/2 fractional quantum Hall 
state 11 - 18 ' and some engineered TSCs such as the surface of 
a topological insulator with an s-wave pair potential 19 ' and a 
semiconductor heterostructure with an s-wave SC and a fer- 
romagnet. 20 ' The latter engineered TSC, which is described 
as a two-dimensional electron gas (2DEG) with Rashba-type 
spin orbit coupling (SOC), Zeeman coupling, and an s-wave 
pair potential and we call the topological s-wave SC, has rota¬ 
tional symmetry derived from the 2DEG realized in the semi¬ 
conductor heterostructure. We expect that the impurity effects 
dependent on the symmetry will be observed more clearly in 
this system. 

There are several ways of treating the impurity effects on 
localized states in a vortex core. The impurity strength T n 
should be small such that the low energy spectra are dis¬ 
crete, since our motivation is to discuss the level spacing in 
the presence of impurities. In this case, we cannot use the 
quasiclassical theory, in which the spectra are treated as con¬ 
tinuous. 21 - 22 ' Moreover, the Born parameter should be van- 
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ishingly small. In this case, there are many but weak scat¬ 
tered, and the (self-consistent) Born approximation is valid. 
In a typical unconventional superconductor, the density of 
states (DoS) is formed around the Fermi energy if the Born 
parameter is not small. 23 ’ Even in an s-wave SC, the impu¬ 
rities, the Born parameter of which is not small, lead to the 
level mixing and the Landau-Zener transition. 24-261 Consider¬ 
ing the above discussion, we can calculate the impurity effects 
using the Kopnin-Kravtsov scheme, 27 ’ in which Green’s func¬ 
tion is used for the CdGM mode with impurity self energy 
[self-consistent Born approximation (SCBA)] while keeping 
the levels discrete. We have improved their scheme in terms of 
the coherence factor and applicability to various types of SC 
in addition to an s-wave SC and call this scheme the improved 
Kopnin-Kravtsov (iKK) scheme. 28 ’ 

The aim of this work is to understand the impurity effects 
on the level spacing related to the adiabatic condition and the 
physical picture of these impurity effects. To address these is¬ 
sues, we first numerically calculate the scattering rates of the 
vortex core states based on the iKK scheme and evaluate the 
property of the minigap for a topological s-wave SC in the 
presence of impurities, in terms of two inequivalent vortices. 
In this case, the sign of chirality related to the chiral edge 
mode is determined by the sign of the Zeeman coupling. 29 ’ 
We find that the obtained numerical results are more compli¬ 
cated than those for the chiral p-wave SC, which have been 
understood only by considering the type of vortex. We then 
make use of the Andreev approximation’ 0-32 ’ in order to un¬ 
derstand the physical picture of these results and the origin of 
the scattering rates. The physical picture can be successfully 
understood by considering the combinations of two types of 
vortex for the chiral p-wave SC. Although the Andreev ap¬ 
proximation is a type of quasiclassical theory, we can use it to 
obtain the physical picture after numerical calculation based 
on the quantum theory. We take ti = 1 throughout this paper. 


2. Model and Method 

We consider the two-dimensional superconducting system 
described by the following Hamiltonian: 19 ’ 

H = J' dr\f{r)H 0 4i{r) + J dr [A(r)^(r)i/rJ(r) + h.c.], (1) 


H 0 = 


P_ 

2m 


- p + a(p x &) z + V z &, 


, <A O) = (t/c|(r),</rj(r)), 

( 2 ) 


where p, a, and V- denote the chemical potential, the strength 
of Rashba SOC, and the Zeeman coupling, respectively. & x ,y, z 
are 2 by 2 Pauli matrices and the symbol denotes a two- 
component vector. The second term on the right-hand side of 
Eq. (1) describes s-wave superconductivity. We consider an 
isolated vortex in the center of a two-dimensional disc with 
radius R and assume A(r) = A(r)e IK8 , where k is the vorticity 
and takes the value ±1. We also assume that A(r) approaches 
a constant value Aq(> 0) far from the vortex core. When V z > 
p 2 + A 2 is satisfied, the system is in the topological phase 
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Fig. 1 . (Color online) Relationships between spectra, impurity scattering 
rates, and minigap are shown schematically. 


and Majorana zero-energy states appear in the vortex core and 
near the edge. From the above Hamiltonian, we can obtain the 
following BdG equation: 

H B d G (r)u K (r ) = E K u K (r), (3) 


#BdG( r ) - 


Ho 
A Hr) 



A (r) = A (r)ifr y . 


(4) 


Here, uk = x {uk‘\,uk\,,Vk j \,Vk\) is a four-component vec¬ 
tor. The angular momenta are good quantum numbers since 
the system has rotational symmetry. The eigenvector can be 
decomposed into angular and radial parts as u K= (j v) (r) = 
Ui(d)u Lv (r) with U,{6) = diag(e®, e i(l+l)e , e ae , e'^ 1 ’ 9 )/ V2^ 
for k = 1 and f/,(6>) = diag(e !(7 " 1 ’«, e ,w , e ,(l+V)e , e ,w ) / V^r for 
k — — 1. Here, we take v as the radial quantum number. We can 
expand this radial part of the eigenstate w/, v (r) by the Fourier- 
Bessel expansion 6 ’ and diagonalize the matrix for each angu¬ 
lar momentum l to obtain the sets of eigenvalues and eigen¬ 
vectors. (Subscripts of the matrix represent the indices of the 
zeroth points of Bessel functions.) We find that there are two 
low-energy modes with energies below the gap in the bulk 
A),; one is localized in the vortex core and the other is near the 
edge. We label them as v = c and v = e respectively. (The def¬ 
inition of Ab is mentioned in Sect. 3.) We remark on the two 
zero-energy states. In this numerical calculation, we consider 
finite-size but sufficiently large systems; hence, the two zero- 
energy states have finite but small energies and their signs are 
opposite. We specify them by v = + and v = —. Each wave 
function has localized distributions around the vortex core and 
the edge simultaneously. We can divide them by taking a suit¬ 
able linear combination of these two states and call the state 
bound in the vortex core v — c. 

We calculate the scattering rates of the vortex core states 
within the single-mode approximation using only the mode 
v = c obtained above. We consider that there exist weak but 
many scatterers and treat them within the SCBA. The scheme 
used to calculate the scattering rates is given in the following 
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form as discussed in Ref. 28: 


c;,r - 1 = Xtt 
r " 2 


(TlticJn) = 


:M/.c(r)(M/, c (r / )) t 
El, c - CTl(iOJ„) - itJn ’ 

M u , 


2n 2 N a ^ E v , c - cri'(ico n ) - ia> n ’ 
M W = j" rdr\(iti, c (r)) f f z iti-, c (r)\ , 


(5) 

( 6 ) 
(7) 


where f, = diag(l, 1, —1, —1) and ui, c (r) denotes the radial 
part of u/, c (r). F n and /V n are, respectively, the impurity scat¬ 
tering rates and the DoS per spin at the Fermi energy ep 
in the normal state. u> n is the fermion Matsubara frequency, 
and in the following, we perform the analytical continuation 
i(o n —» (i) + id. The DoS can be obtained from Green’s func¬ 
tion: 


l 


drlm —Tr f z G(r, r, ioj n )| 


icon — 


( 8 ) 


We define the impurity scattering rates (denoted as T) in two 
ways. The first definition of T is the half width at half maxi¬ 
mum of the spectrum, shown as F in Fig. 1, and we use this 
definition in Sect. 3. The second definition of T is the inverse 
of the DoS multiplied by 7r [T = (ImG)" -1 ], shown as T' in 
Fig. 1, and we use this definition in Sect. 4. These definitions 
are equivalent when the spectrum takes the Lorentzian profile. 


3. Numerical Calculation 

In this section, we show the results of numerical calcu¬ 
lations for scattering rates. First, we explain the parameters 
used in our calculations. We fix the quasiclassical parameter 
= 20, where kt and £' are the Fermi momentum and co¬ 
herence length, respectively. In order to set this parameter, we 
consider the homogeneous system with the s-wave pair poten¬ 
tial A (r) = Ao for the moment. In the topological phase, the 
Fermi momentum kt and the Fermi velocity vp are well de¬ 
fined since this two-band model has a single Fermi surface for 
the uniform system. We set Ab as the minimum gap of the one- 
particle excitation spectrum around k = kp and define the co¬ 
herence length as £ = vp/ Ab with this minimum excitation gap 
Ab. [Note the difference between Ao and Ab; Ay is simply the 
amplitude of the s-wave pair potential and Ab is the minimum 
excitation gap. In general, they are not the same and approxi¬ 
mately satisfy |Ab/Ay| - \a\k?l(V} + a 2 kp )^ 2 as Eq. (34).] We 
set ii = 0 and Ao as the unit of energy and change the ratio of 
ma 2 /\V z \ under the condition V 2 > /j 2 + A 2 . We remark again 
that this system has two inequivalent vortices in terms of the 
symmetry operation. We can distinguish this inequivalence by 
the relative sign of the slopes of the vortex core mode and chi¬ 
ral edge mode: dE[,Jdl and dE/, e /dl, as shown in Fig. 2. We 
also find that the signs of the Zeeman coupling V, and the 
vorticity correspond to the signs of the slopes of the chiral 
edge mode and vortex core mode, respectively. 29 ’ This struc¬ 
ture is very similar to that of the chiral p-wave SC if we regard 
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Fig. 2. (Color online) Energy spectrum for fi = 0, = 20, and 

ma 2 /\V z \ = 1.6. (a) The sign of the Zeeman coupling is positive and the 
slopes of the two low-energy branches are opposite in sign, (b) The sign of 
the Zeeman coupling is negative and the slopes have the same sign. The blue 
solid line labeled as “Andreev” is the analytical solution for the CdGM mode 
discussed in Sect. 4. 


the sign of the slope of the edge mode as the chirality of the 
Cooper pair. We thus call the type of vortex with the positive 
(negative) relative sign the “parallel (antiparallel) vortex”. In 
this system, the signs of the two slopes are the same (oppo¬ 
site) when kV z > 0 (kV 7 < 0), and in the following numerical 
calculation, we only change the sign of the Zeeman coupling 
and fix the vorticity as k — -1. To solve the BdG equation, we 
set the system size to R = 20k and the spatial profile of the 
pak potential in the radial direction to A(r) = Ay tanh(r/£). 
Moreover, we introduce two cutoffs, l c = 50 and N c = 400, 
which are the numbers of angular momentum as a quantum 
number and zero points of the Bessel function, respectively, 
and the infinitesimal quantity 6 = 10 6 A h in Eq. (8). Under 
these parameters, we numerically solve the BdG equation and 
Eqs. (5) - (8) to calculate the DoS and scattering rates in the 
presence of impurities. 

We show the scattering rates T for various values of 
ma 2 /\V z \ in Figs. 3(a) and 3(b), which correspond to the pos¬ 
itive and negative Zeeman couplings, respectively. The hori¬ 
zontal axis shows the eigenenergy for pure systems, which is 
scaled by the gap Ab for each parameter. The discrete eigenen- 
ergies scaled by Ab are independent of the ratio ma 2 /V z , and 
the level spacings scaled by Ab in the low-energy region are 
the inverse of k^. We set T n = l(T 3 ;rAb. Note again that Ab is 
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Fig. 3. (Color online) Impurity scattering rates for various values of 
ma 2 /V z and T n /7rAb = 10~ 3 . (a) V z > 0 (antiparallel vortex) and (b) V z < 0 
(parallel vortex). 


determined by some parameters such as V z , ma 2 , and so on. 

These figures show that the scattering rates are smaller for 
V z > 0 than for V z < 0. One might consider that the impurity 
effects are characterized by the type of vortex (parallel or an¬ 
tiparallel), because in the chiral p-wave SC, the bound states 
in an antiparallel vortex are robust against impurities, while 
those in a parallel vortex are sensitive to impurities. It is im¬ 
possible to characterize the impurity effects only by means of 
the relative sign of the vorticity and chirality since there are fi¬ 
nite scattering rates even for the antiparallel vortex, as shown 
in Fig. 3(a). We need to consider the dependence of scattering 
rates on another parameter ma 2 /V z . We find that the scattering 
rates in the low-energy region are a decreasing function of the 
ratio ma 2 /\V Z \ for the parallel vortex [Fig. 3(b)], while they 
are increasing function for the antiparallel vortex [Fig. 3(a)]. 
In particular, for very small ma 2 /\V z \ of 0.05 and 0.1, the scat¬ 
tering rates are exceptionally suppressed for the antiparallel 
vortex, and the difference between the two inequivalent vor¬ 
tices is evident, similarly to the chiral p-wave SC. On the other 
hand, when ma 2 /\V z \ is large, the scattering rates of the two 
vortices have almost the same energy dependence and mag¬ 
nitude, which is a different feature from the case of the chiral 
p-wave SC. (We comment on the exceptionally small scatter¬ 
ing rates of zero-energy states in spite of the finite scattering 
rates of other excited states. This is independent of the type of 
vortex or the parameter ma 2 /V z .) 



Fig. 4. (Color online) Impurity effects on minigap. The horizontal axis de¬ 
notes the impurity strength r n /(/rAb). The symbols have the same meanings 
as those in Fig. 3. 


We note again that the adiabaticity of TQC requires the ro¬ 
bustness of minigaps, and thus we show the F n dependence 
of the minigap in Fig. 4(a) for V z > 0 and in Fig. 4(b) for 
V z < 0. We define the minigap in the presence of impuri¬ 
ties by subtracting the widths of the first excited state and 
zero mode from the minigap for a pure system, as shown in 
Fig. 1. Impurity effects on minigaps reflect the scattering rates 
in low-energy states. In Fig. 4(a), only for the case of small 
ma 2 /\V z \ with V z > 0, the minigap is robust against the im¬ 
purities; otherwise, the size of the minigap steeply decreases 
with increasing F n . In this system, the stability of the minigap 
crucially depends on the type of vortices and the parameter 
ma 2 /V z . 

We consider the origin of the dependence of the scattering 
rates on the type of vortex and the parameter within the An¬ 
dreev approximation in the next section. 

4. Andreev Approximation 

4.1 Spectrum and wave functions of low-energy states 

In this section, we analytically study the bound states in the 
vortex core within the Andreev approximation for the BdG 
equation. Tewari et al. also discussed the zero-energy states 
in the vortex core by solving the approximated BdG equa¬ 
tion analytically. 32) The main strategy of the approximation 
is the same as in the following discussion, but the details of 
the calculation are slightly different. Moreover, our aim is not 
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only the construction of a zero-energy state but also the con¬ 
struction of nonzero-energy states and the understanding of 
the scattering rates. In Appendix A, we also discuss the possi¬ 
bility of the construction of the edge state within the Andreev 
approximation. 

We start with the BdG equation for the uniform case. 


HBdG(k)i<k = Euk , 

where Hsddk) is a 4 by 4 matrix and described as 

A 


HBdG(k) = 


H 0 (k) 
A 4 


-H*(-k) 


H 0 (k) = 


ft- + v z 
-iake'^ k 


iake~ ilpk 
e k - V z 


A = 


0 

-A 


H = - - R, 

2m 


k = kcosi(f>k)e x + ksm(<pk)e y 


(9) 


( 10 ) 


( 11 ) 


( 12 ) 


ilk is a four-component vector and taken as u k = e lk r Uk . Uk is 
a spatially oscillating wave function. 

First of all, we take the unitary transformation ii k = 
U(k)it k = e lkr Uik)u t j c , in which basis Hoik) is a diagonal 
matrix. One of the unitary matrices is taken as 

Wiik) 6 


U(k) = 


Ufk) = 


[ 

6 u\i-k) \ ’ 


1 

iake -i *‘ 

J\H + a 2 k 2 - V z 

c 

Jv* + a 2 k 2 - V z 

iake'* k 


(13) 


(14) 


C~ = 2 


^V} + a 2 k 2 v (+ a 2 k 2 v - l/ z ). (15) 

Through this transformation, we obtain the following Hamil¬ 
tonian: 33 ^ 


Hi, 


J BdG (k)u k - Eu k , 
H h BdG ik) = uUk)H Bd Gik)U(k) = 

diag iE k+ ,E k J),E k± = e k + ^V 2 + a 2 k 2 , 



A b ‘ 


-H b 0 _ 


% = < 


A b = 


fsik) 


V, 


-ifpik) Ae‘* k f s ik) A 

-Uk) A if p ik) Ae~» k 

ak 


f P ik) 


Jvl + a 2 k 2 Jv - + a 2 k 2 


(16) 

(17) 

(18) 

(19) 

( 20 ) 


E k+ - E kf+ = 2 Jv* + a 2 k 2 F = E c , 


( 21 ) 


,, dE k . 

Ek- - -ITT- ■ (k - k F ) 
ak 

2 

mcv 

i - 


---- ik - k F ) = v F -(k- k F ). 

1 »■ 


( 22 ) 

Here, the Fermi velocity v F is parallel to the Fermi momentum 
kp. In the following part, we consider a single vortex in the 
system and set A = A (r)e ,Ke (/r: vorticity). We define a 2 by 
2 matrix Rid) as the rotation around the z axis by angle 6 and 
introduce the basis of the 2D polar coordinates as (e r , eg) = 
R(0)ie. x ,e y ). In the bulk region far from the vortex core, A(r) 
approaches Ao(> 0), which is the magnitude of the s-wave 
pair potential. Because of this inhomogeneity, we replace (k- 
k F ) by —/V in Eq. (22). In the following, we simply use f p and 
f s instead of f p ikp) and f s ikp), respectively. We can thus write 
down the differential equation for the slowly varying function 
u° k as 

Ecu b h+ - ifpe** Av b p+ + /,Av*. _ = Eu b k¥+ , (23) 
-/v F • VS b p _ - /;Av b p+ + ifpe-V* Av b kF _ = Eu\ ¥ _, (24) 

-Ecv* + ifpe-^X 


< + 


• f.A*u h 


kp— - E^k ¥+ , 


(25) 


fv F • Vv£ p _ + / s A‘« b p+ - */ P ^ATl b p _ = £v b p _. (26) 

We can easily eliminate u° kf+ and v^ + and obtain the follow¬ 
ing two coupled differential equations: 

-iv p • Vcr z - ReA p 2 <T v + lrciA pl & y 




+A,AA s j |pb F- j - E\ 
where we define A p i 2 = -if p e ±u,>k r A, A s = f s A and 


u°. 


,b F ~ 

kp~ 


(27) 


A = 


e 2 - El - |A p ii 2 [ a;, 


E + E r 


A p i 

E-Er 


A, = —i 


0 

-A! 


A s 

0 


For V 2 > A 2 + /r 2 , we find E k+ > 0 and the existence of 
kp such that E kv - — 0. We expand the above BdG equation 
[Eq. (16)] around kp while retaining the leading-order terms 
only, i.e., up to the first order for E k - and the zeroth order for 
the other elements. We call this approximation the Andreev 
approximation. We introduce two quantities as follows: 


We can construct the zero-energy state (E = 0) analogous to 
the Jackiw-Rebbi solution for the massive Dirac equation 34 ' 
and obtain the low-energy spectrum through the first-order 
perturbation even in the presence of A. However, the contri¬ 
bution from A, i.e., (M b p+ , v b p+ ), is small compared with that 
from (n b p , v b p _) and we can thus omit the fourth term on the 
left-hand side, of Eq. (27) and hereafter analyze the following 
equation: 

[-'V F • Vd-, - ReA p2 <x v + ImA p2 r> v | = El fe F 

\ kp— / \ kp— 

(28) 

We introduce another set of coordinates (,v, b) as (e s . eg) — 
R((l>k)(e x , <?v). In this frame, Vp = vpe s . We take 
the gauge transformation as *(M b 


r,b 


U f -> V *p-) = eX P['^ 
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l)0t F (T,/2] t (Mj. F _, v£ p _); then Eq. (28) is reduced to 


-iv F d s fr z - f p Mr)-&y - K f p A(r)-cr_ 


y *Jp l 
r r 



(29) 


Note that r = s/s 1 + b 1 and tan (8-(f>k F ) - b/s. b is the impact 
parameter in scattering theory and related to the angular mo¬ 
mentum l through / = -k F b. From this point, we find that there 
exists a zero-energy state when b = 0. For small b , we treat 
the third term of Eq. (29) as a perturbation Hamiltonian and 
obtain other low-energy states. We can construct the solution 
with the energy E - 0 for the unperturbed Hamiltonian: 

| : ) = exp[- ! « s .l,)](_ sg 1 ii(ii) ). (30) 



where r' = Vs ' 2 + b 2 . With this wave function, we can obtain 
the correction of energy through the first-order perturbation: 

/-*oo 

E(b) = AEi = I dsK\f p \ A(r)(b/r) exp [-2 ip(s, b )] /Nf, 
Jo 

(32) 


N 


=-f 


ds exp [- 2 /?)] /£. 


(33) 


The solid lines shown in Fig. 2 represent the energy spec¬ 
trum given by Eq. (32) for kyj = 20 and ma : /\V z \ = 1.6. For 
this parameter, Eq. (32) is in good agreement with the energy 
spectrum numerically obtained from the BdG equation. For a 
smaller mcr/ |V.|, the deviation of the energy spectrum for a 
large angular momentum is larger, but the low-energy spec¬ 
trum can be described well. From this, an expression for the 
minigap can be approximately obtained: 


Amini — 


1 dE 


k F db 

b =0 


Cl¬ 


ip Ao 


Cl 


Ah 

k F % 


(34) 


where ci is a constant on the order of unity and, in particular, 
ci = 7£(3)/n 2 when A (r) = Aotanh(r/£). 


4.2 Scattering rates 

Hereafter, making use of the above solutions to the Andreev 
equation, we calculate the scattering rates for low-energy 
states. In the following calculation, as mentioned in Sect. 2, 
we change the definition of T to the inverse of the DoS multi¬ 
plied by n [T = (ImG)^ 1 ]. We start with Eq. ( 6 ). We note that 
the magnitude of the momentum k is k F within the Andreev 
approximation. Hence, we omit the subscript of </i/ ;| and also 
use tp instead of k F as a label. We define ip , which satisfies 
ip = 26 - cp + n\ s = -s and b = b. Without loss of generality, 
we can take </> such that \9 - <p\ < n/2 and, in this case, s > 0. 
In order to use Eq. ( 6 ), it is necessary to construct an approx¬ 
imate eigenfunction for the BdG equation. Here, we assume 
that we can approximately describe the eigenfunction using 


the following formula: 

e iW u,p(s, b ) + c 2 Uj,(s, b) 

uiJr) = —p -r---, (35) 

yj2nNd^r 

where d is the number of Fermi surfaces in the normal state 
and c 2 is the phase factor discussed in Appendix B. We can 
confirm that this formula is valid for an s-wave SC and a 
chiral p-wave SC through their analytical expression for the 
BdG wave function (we directly confirm this in Appendix B). 
We need to obtain u$ from through the unitary and gauge 
transformations since 14 used in Eq. (35) is described in the 
original basis. Moreover, we assume that the matrix element 
s u^T z u<p 2 satisfies the relation = —m ( p 2 _ l p l . This 

relation is satisfied in the present system, an s-wave SC, and a 
chiral p-wave SC. Under these assumptions, we evaluate the 
diagonal elements of Eq. (7): 


Mu = rdr\lTiu) {r)T z ui c (r)\- 


S 

-f 

■f 




_ e -2% s + c * m ^ g 


2ikps\ 


r (4 Nd& 2 

I |2 

dr \m^\ [1 - cos(4fc F s - 2/3)] 
~ 8 (Nd^) 2 ' 


(36) 


In the last line, we put c 2 = e'P, and when k F ^ is large, the 
integral of the second term is negligible. This expression for 
the matrix element described by the Andreev wave function is 
general and we also calculate the scattering rates in the case of 
the s-wave SC using the formula in Appendix B. We evaluate 
the matrix element Im^l 2 as 



jv?+a 2 k 2 F -V z ^V 2 +a 2 k 2 F +V z 

y jv 2 + a 2 k F ^jv 2 + a 2 k F 

x |sin (0 - ip )| exp [- 2 i^(s, b)\. 


S K . 1 
(37) 


This expression shows that the impurity effects are the same 
even if we change the sign of V z and k simultaneously. The 
difference between parallel and antiparallel vortices is de¬ 
scribed only by the coefficients [(V 2 + a 2 k 2 F )^ 2 + |Vj|]/(V? + 
a 2 k 2 F )^ 2 [+ (-): parallel (antiparallel)]. In order to evaluate 
the integral, we use 

2\sb\ 

| sin(0 - (p)\ = 2\ sin(0 - (p ) cos (9 - (p)\ = — ——, (38) 

s~ + b~ 

and perform the following approximation in the region of the 
integral. The integral in Eq. (36) can be evaluated as 



dr . 7 
— sin ~(0 - <p)e 
r 


-4 ijj(s,b) 


dr 4s 2 b 2 

J h r (s 2 + b 2 ) 4 


= 1 - 


2 b 2 b 4 

l r + ¥~ 


(39) 


If we focus on the first excited state with angular momentum 
|/| = 1 , we can estimate b 2 /£ 2 = l/(kF^) 2 and omit the second 
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and third terms in Eq. (39). We thus obtain the matrix element 
Mi l as 


M U = 


1 


8 (Nd£) 2 


1 + 


sgn(ArVj) \V Z \ 




• a 2 kl 


(40) 


We consider the two limiting cases of the scattering rates in 
detail. One is that there is a large level spacing in the pres¬ 
ence of impurities and hence we can neglect the contribution 
from other angular momentum states. The other is the limit 
of continuous spectra, but the impurity strength is very small 
and thus we can treat it within the non-SCBA. 

In the former case, we neglect the contribution to 07 from 
other angular momentum states l' in Eq. ( 6 ); hence. 


cr,(ioj n ) = 


T n _ M U 

2ji 2 N „ E Uc - <Tl(iu ) n ) - iio n ’ 


(41) 


cri(E Uc ) = i 



(42) 


E = Imcr i(Eic) 


= r n 


7rA b 

f \ 

, , sgn(/eK) \V Z \ 

8(7rNd) 2 r n kp^ 

V 

^Vr + a 2 k 2 


(43) 


In the second line, we perform the analytical continuation 
ia>„ —> to + id and set u> = Ej tC . We also use N n = kp/(2nvp) in 
the last line. 

In the latter case, we replace the summation with respect 
to /' by integration with respect to b' through the relation 
l' = -kpb'. We focus on the low-energy states and thus use 
the approximate form of the energy spectrum £;- c = E(b') - 
KA umu k[-b' and put <t>( iu> n ) - 0 in the denominator in Eq. ( 6 ), 
and then we can calculate T as 


r = lmcr,(E(b)) = 


kfT n 

2n 2 N n 


f 


db'lm- 


Mi(b)Kf) 


E(b ') - EQ}) - id 


r n M)j 

2nN n A lrnn j 


r„ 


8ci (Nd) 2 


t ' 

+ sgn(*V.) |W| 

y jv 2 7a 2 l^ / 


(44) 


We compare the above two formulas with the numerical 
calculation. We again comment on the definition of scattering 
rates. In this section, in contrast to the previous section, we 
define T as the inverse of ImG at the energy £) c . In Figs. 5(a) 
and 5(b), we show the scattering rates for the first excited 
state with angular momentum l - 1. The solid lines shown 
in Figs. 5(a) and 5(b) represent the scattering rates calculated 
using Eq. (43) while the dotted lines shown in Fig. 5(b) rep¬ 
resent those calculated using Eq. (44). The dashed lines with 
open squares and open circles in Figs. 5(a) and 5(b) repre¬ 
sent the scattering rates based on the iKK scheme represented 
by Eqs. (5) - (7) for V z > 0 and V z < 0, respectively. All 
the schemes show that when V z > 0, i.e., an antiparallel vor¬ 
tex, the scattering rate is an increasing function of ma 2 /\V z \, 
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ma 2 /\V Z \ 



Fig. 5. (Color online) Scattering rates calculated by means of Andreev ap¬ 
proximation [the solid lines labeled as '‘Andreev( 1)” correspond to the calcu¬ 
lation using Eq. (43) and the dotted lines labeled as “Andreev(2)” correspond 
to that using Eq. (44)] and the scheme represented by Eqs. (5) - (7) (open 
squares for V- > 0 and open circles for V z < 0). 
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Fig. 6. (Color online) r„ dependence of scattering rate 1. The solid line 
denotes the behavior of [r n /(ffAb)E^ 2 . 


and when V z < 0, i.e., a parallel vortex, the scattering rate 
is a decreasing function of ma 2 /\V z \ as mentioned in Sect. 3. 
We need to explain the vertical axis in Figs. 5(a) and 5(b). 
In Fig. 5(a), we set T n /( 7 rAb) = 10 4 . The magnitude of the 
scattering rates based on the iKK scheme is shown on the 
left vertical axis and that based on Eq. (43) is shown on the 
right vertical axis, the range of which is 0.64 times that of 
the left axis. We tune the range in order to make the param- 
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eter dependence clear. In Fig. 5(b), we set r n /(7rA h ) = 10 1 . 
F based on the iKK scheme is shown on the left axis and that 
based on Eq. (43) is shown on the right axis, the range of 
which is 0.33 times that of the left axis, as well as in Fig. 5(a). 
Moreover, I calculated by multiplying Eq. (44) by V2 is also 
shown on the left axis. (The aim of this operation is simply 
to make the parameter dependence clearer.) The level spac¬ 
ing in this case is A m ; n i/Ab ~ (£f(t) 1 = 0.05, which is much 
larger than r/A b ~ 4-n x ICC 4 in Fig. 5(a). In this case, we 
can consider that the system is in the limit of discrete spec¬ 
tra and that Eq. (43) is valid. On the other hand, in Fig. 5(b), 
T/Ab ~ 0.27T x 10 1 , and this is comparable to the level spac¬ 
ing. Hence, Eq. (43) no longer describes the ma 2 /V z depen¬ 
dence of T /r n while this dependence is similar to Eq. (44), ex¬ 
cept in the region of very small ma 2 /\V Z \ for V- > 0. In this ex¬ 
ceptional region, the system has discrete spectra, i.e., the level 
spacing A m ; n i is always larger than the scattering rate F since 
T —» 0 as ma 2 /V z —> 0. We infer that Eq. (44) does not explain 
the dependence precisely even when V z < 0 for the following 
two reasons. First, we cannot regard the system as being in 
the continuous limit because the scattering rate I is not much 
larger than and comparable to the level spacing A m ; n i. Second, 
we should treat the impurities within the SCBA rather than 
the non-SCBA because the impurity strength I n is large com¬ 
pared with the level spacing. We emphasize, however, that the 
ma 2 1 IV-1 dependence of F/r n is similar to Eq. (44) rather than 
Eq. (43). 

Figure 6 shows the F n /(7rAb) dependence of r/r n . We ex¬ 
pect r /r n oc 1 / Vr n /(7rA b ) from Eq. (43) and Fig. 6 implies 
that this is a good description for small r n /(^A b ), approxi¬ 
mately for r n /(7rA b ) such that < 0.1. 

From the above discussion, we find that it is possible to 
describe impurity effects by the Andreev approximation. In 
the remaining part of this article, we discuss the origin of the 
scattering rates in terms of the “spin-resolved angular mo¬ 
mentum” defined below and connect it to the “total angular 
momentum”, which is the sum of the vorticity and chirality. 
As mentioned in the text just below Eq. (4) in Sect. 2, each 
component of the quasiparticle state labeled as l has the angu¬ 
lar momenta (/-1 T, Z /+1 T, /1) and (/ 1,1+1 |,/T,/-1 I) 
in the original basis when k — -1 and 1, respectively. We 
can understand the origin of the scattering rates by dividing 
these four components into two pairs such that ( ler , ler) and 
(l + ler', l + ler'). We introduce the spin-resolved angular mo¬ 
mentum to label these pairs: 

L zcr = -26 K -i6a-,-i + 2 < 5 k i i < 5 0 - i |. ( 45 ) 

Here, we review the calculation of Eq. (37) in terms of L zcr . 
We symbolically describe a part of the unitary transformation 
u —> u° as 

u b - - ^ a-o-Wo-, v b _ = ^ b-o-Vo- (46) 

<X (T 

We note that \a-o\ = \b-a\ and |a_-f| 2 + |a_jj 2 = 1. The matrix 


element can also be described as 

= |a~oj 2 |sin (0 - (j>)\. 

(47) 

cr denotes the spin component with \L z g\ = 2 while the other 
spin component with L zcr = 0 does not contribute to the scat¬ 
tering rates. We remark that | a-o-\(cr =T, I) are the degrees of 
mixing of iq and tq to make the Fermi surface. Therefore, the 
scattering rates are determined by the ratio of the contribution 
from the spin component with \L za \ = 2 to the Fermi surface. 

Finally, we introduce the total angular momentum L z = 
-i2f z (dg+d,p), which acts on the wave function e~ ll6 uijr). We 
can corroborate that this definition is consistent with Eq. (2.9) 
in Ref. 35 for the s-wave SC and chiral p-wave SC. We calcu¬ 
late the expectation value 

(L z )=Tr\til c (r)e il %e- i %Ar)\ 

= Y \a-o\ 2 L Z o- = \a- & \ 2 L z& . (48) 

cr 

Therefore, we can understand the impurity scattering rates as 
the magnitude of the total angular momentum. Both limiting 
cases, | (L z ) \ = 0 and 2, where L- is a good quantum num¬ 
ber, correspond to chiral p-wave SCs. On the other hand, it is 
remarkable that | (L z ) \ = 1 does not describe the s-wave SC 
because L z is no longer a good quantum number. | (L z ) \ = 1 is 
simply due to the superposition of | (L z ) \ = 0 and [ (L z ) \ = 2 

5. Summary 

In this paper, we have investigated the impurity effects on 
the bound (CdGM) states in a vortex core of a topological s- 
wave SC. We have calculated the scattering rates for CdGM 
states numerically and analytically. The numerical calculation 
is performed for discrete spectra based on the iKK scheme in 
order to discuss the impurity effects on the level spacing of 
bound states. We have found that there are similarities and 
differences compared with the chiral p-wave SC: the parame¬ 
ter ma 2 /V z is important for describing the impurity effects as 
well as the relative sign of the vorticity and chirality. On the 
other hand, the analytical calculation is based on the Andreev 
approximation to understand the physical picture of impurity 
scattering and the origin of the parameter-dependent scatter¬ 
ing rates. The conclusion is that the scattering rates \a-o\ 2 are 
determined by the ratio of the spin components with \L za \ = 2 
composing the Fermi surface, and this ratio describes the 
magnitude of the total angular momentum | (L z ) \ = 2\a- cr \ 2 . 

Although the topological s-wave SC belongs to the same 
class as the chiral p-wave SC, the properties of excited states 
are more complicated than those in the chiral p-wave SC. 
However, vanishingly small scattering rates because of the 
coherence factors and rotational symmetry are obtained as 
well as in the case of the chiral p-wave SC. The topological 
s-wave SC is more promising than a chiral p-wave SC such 
as Sr 2 Ru 04 in terms of the robust low-energy states in a vor¬ 
tex core in TSC owing to the continuous rotational symme- 


^|a_ ;r | 2 sin^^(0-0)j 
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try since this superconductivity is expected to appear in the 
2DEG in the semiconductor heterostructure while the chiral 
p-wave superconductivity in Sr 2 Ru 04 is affected by the C 4 
point group symmetry. 
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values and eigenvectors: 

aM = ^V}-/u 2 -Aji= \-b + , -b., b + , bJ), (A- 8 ) 

aA 2 = yjvd-i-i 2 - A, b 2 = \b + , b-, b + , b.), (A-9) 

b + = y/V~TJi/2\V z \, b- = JV^Ti/ 2\V z \. (A-10) 

With these sets of eigenvalues and eigenvectors, wave func¬ 
tions can be constructed as 


Appendix A: Edge state 

We construct the zero-energy edge state, which appears 
near the edge of a TSC, in order to confirm that the Andreev 
approximation can describe the features of a TSC. For sim¬ 
plicity, we consider the one-dimensional system of a topolog¬ 
ical s-wave SC in the region x > 0 and vacuum in the region 
x < 0. The BdG equation for the homogeneous topological 
s-wave SC in one dimension is as follows: 


H(k)u k = Eu k , (A-l) 

H(k ) = ekf z + V z & z f z - ak& y f z - A& y f y , (A-2) 

where k is the one-dimensional momentum. & and f denote 
Pauli matrices in spin space and Nambu space, respectively, 
and they are 4 by 4 matrices. The vector u k has four com¬ 
ponents, Uk = L (uk], Uki, Vfcf, In a similar way to Sect. 4, 
we can write down the Andreev equation by introducing a 
slow spatial variance as an edge at x = 0 , but in this case, 
we impose the boundary condition where the wave function it 
is 0. Our goal is to construct the low-energy, particularly the 
zero-energy wave function satisfying the boundary condition 
by linear combination of the zero-energy solutions to the An¬ 
dreev equation. The Andreev equation with k = ±k F can be 
written as 


+ | md x (f z + f p A(x)& y \ u° ±kf _ = 0. 
The solutions to these equations are given by 
1 


(A-3) 


~/b 

U 


.±kc 


exp 


1 r x 

ikpx - f p A(x')dx' 

Vf Jo 


(A-4) 


The unitary matrix Lf defined as 14 = UkU° k is given by 


U k = 


(V 


iakl + JV 2 + a 2 k 2 - V, 


/C, 


(A-5) 


where C is defined by Eq. (15). In addition to these solutions 
with Fermi momenta, we construct the zero-energy solutions 
with k ~ 0. In this case, we can write down the Andreev equa¬ 
tion as 


\V z & z t z - /jt z + iad x &yf z - Ao-yf^J uq = 0. (A- 6 ) 

This eigenvectors can be obtained by considering the matrix 

B = \ V z & x + ifi(f y + Af c ] /a. (A-7) 

We can construct the zero-energy solutions with positive 
eigenvalues A,(i = 1,2, • • •) and the corresponding eigenvec¬ 
tors b, of B. When V 2 > A 2 + /.r, there are two positive eigen- 


it 0 /i, = bi exp 



(i = 1,2), 


(A-11) 


which decay far from the edge at x = 0. We consider the linear 
combination of these bases: 


it(x) = C + u + k ¥ + C-it-ky + CiitoA, + C 2 uoa 2 - (A- 12 ) 

Here, we omit the contribution from it ±kv+ and perform the 
unitary transformation U ± k r to obtain u ± k r . Since we can 
choose C+ and C such that C+it+ky + C-U-ky is parallel to 
i7( U at x = 0, we find that C 2 = 0. We emphasize that it is 
necessary to consider the k = 0 components only for the edge 
state rather than the vortex core state. 


Appendix B: Scattering rates for s-wave SC 

The BdG wave function for the CdGM states in a vortex 
core for the s-wave or chiral p-wave SC and its asymptotic 
behavior for large kpr are described as 

1 ( Ji-L z ,2(kpr)e- K ^e^- L ^ e 

v ; /,cj.( r )/ Viv 

o il0—K(f) ( --M/2 I -T-ndc-ari . ~—iQi-T-n(kar\ I \ 


Hu = 


g-iS/2 e igi-L-j2(kyr) + e -igi-4/ 2 (fer)] 


yj2nNk F r {-ie iL ^ /2+ie/2 + e -igi^, 2 (kyr) jj 

(B-l) 


K(r) = 


vf 


f 


A (r')dr'. 


(B-2) 


gi(x) = x + l /(2x) - (21 + l)zr/4, (B-3) 

where N is a normalization constant. -L z is the sum of the 
vorticity and chirality and for example, L z takes +1 for the s- 
wave SC, +2 for the chiral p-wave SC with the parallel vortex, 
and 0 for the chiral p-wave SC with the antiparallel vortex. For 
convenience, we adopt the minus sign in the definition of L z . 
Equation (B-l) is reduced to 

/ tq,cf(T)\ 

\V/, 4 (r)/ 


c^e 


y/2nNk F r 


'i[gi,i(kyr)—g 2 ,i(kyr)] 1 ' 
le is l 2 {gUguM+guikyr)] + i c ^ e -i\gu(kyr)+g 2 J (kyrj\^ 


I g-M! 1 e ilgu(kyr)-g 2 j(kyr)] + - c £ 


g u (x) = x + (/ 2 + Z|/ 4)/(2jc), 
gi,i(x) = L z l/( 2x), 


C3 = g -i( 2;-4 + l>/4^ 


, (B-4) 

(B-5) 

(B- 6 ) 

(B-7) 
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C 4 = e 


i(2l-L z )n/2 


(B- 8 ) 


C 3 is simply an overall phase factor and we need not consider 
it. We take an approximation that gi,i(kpr) ~ /cf|,v|. Of course, 
this approximation is not valid in the region where £f|s| is not 
large compared with |/| = kp\b\, but we expect that the physical 
picture does not change significantly. 

On the other hand, the Andreev wave function is described 
as 

(^) = e iL ^/2 j exp [ik F s - 1 Ks, b)], (B-9) 

1 r s s' 

if/(s, b) - — A (r')-ds'. (B-10) 

vf Jo r' 

As mentioned in the main text, we take <p such that \6 - <p\ < 
n1 2, hence, s > 0. We consider a linear combination of two 
momenta <p and <A, which share the impact parameter b, 

e^ b % + c 2 ^\ 

( [gikps+iL z b/(2r) _|_ ^^-iL z jr/2^-ikps—iL z b/(2r)^~iL z 6/2\ 

_i[ e ikw-iLzbl(2r) + c ^ e iL z n/2 e -ik P s+iL z b/(2r)-^ e iL z 6/2 I 

( ^ikpS-iL z l/(2kpr) _j_ ^^-iL z n/2^-ikps+iL z l/(2kpr)^-iL z G/2\ 

_i\gikps+iL z l/(2kpr) _|_ ^^iL z n /2 g—ikps-iL z l/(2kpr)^iL z G /2 J > 

(Bll) 


where we use 6 - <p ^ sin( 6 > - (f>) - b/r in the second line and 
l = -k\-b in the third line. When b is small compared with s, 
K(r ) ~ i//(s, b ), and thus we obtain 


C2 = W4 — e 


i(/+l/2)jr 


(B-12) 


JIB 


—l,c 


sjlnNk F r 


l lf, + C2ip] ■ 


(B-13) 


We can choose the normalization constant N = ANJdIky, 
where N is defined by Eq. (33) and d is the number of Fermi 
surfaces; in this case, d = 1. Equation (35) can thus be ob¬ 
tained in the case of s-wave and chiral p-wave SCs. 

In the remaining part of this appendix, using Eqs. (36) and 
(44), we calculate the scattering rates of the CdGM mode in 
a vortex core in an s-wave SC. In this case, / is a half-odd- 
integer and L- = +1. We evaluate the matrix element Mu 
using Eq. (36), 


M„ = 


-l 

-I 

-JT 


dr \ m !f>ii>\ 

~ m 2 f- 

v dr 12; sin [L z (cf> - 0)/2]|“ 


-4 i/s(s,b) 


h r m 2 f- 

^ scos 2 (6-(f>) 

o dS ~- 2 N 2 f 
3 1 


ds 


0 r 4 2 N 2 ? 


1 


2 N 2 f 


ln \b\ + — 


1 - ln(^ z + b~) ■ 


if- + b 2 


1 


■ln(| b\/0 (as b —> 0 ). 


(B-14) 


2N 2 f 

Noting that \E(b)\ = Aminikpl^l- the scattering rate is obtained 
as 


r 

t: 


1 


47rA 2 ^ 2 A n A n 


■In 


E(b) 


Amini^F^ 


c\N 2 


In 


E(b) 


Amini^F^ 


(B-15) 


where the coefficient is on the order of 1. This logarithmic 
dependence of E on E(b) is in good agreement with earlier 
works based on the quasiclassical theory. 36,37) Similarly, we 
can calculate the scattering rates of CdGM states for chiral 
p-wave vortices. 
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